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Abstract 

We present a calculation of the change of free energy of a solid surface upon 
bending of the solid. It is based on extracting the surface stress through a molecular 
dynamics simulation of a bent slab by using a generalized stress theorem formula, 
and subsequent integration of the stress with respect to strain as a function of 
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bending curvature. The method is exemphfied by obtaining and comparing free 
energy changes with curvature of various reconstructed Au(OOl) surfaces. 

KEYWORDS: Molecular dynamics, Bending of surfaces, Surface thermody- 
namics. Gold, Curved Surfaces. 



1 Introduction 

The surface free energy of a solid, defined as the excess free energy of the semi- 
infinite solid relative to an infinite solid with equal number of particles is 
a very important but also a very elusive quantity. It is important because 
its minimum determines the surface equilibrium state and its properties. 
However, while in the liquid it coincides with the ordinary surface tension, in 
the solid it becomes elusive because there it can neither be easily measured 
nor calculated. The quantity which is easy to access is instead the surface 
stress (generally a rank 2 tensor) defined as the change of the surface free 
energy with respect to a unitary increase of surface area upon stretching, 
which combines surface specific free energy 7 and its first derivative: 

dA 

While 7 is of course positive for a substance below its critical point, the 
magnitude of the second term (arising from the rigidity of the solid) is gener- 
ally comparable, but can have either sign. Therefore, while surface stress 
is becoming increasingly available through measurement , simulation^], 
and microscopic calculation |Q, the surface free energy remains generally un- 



known. At T = 0, where surface free energy and surface energy coincide, 
there are of course a large number of calculations, as quoted e.g. in Ref. 
0). However, finite temperature free energy calculations, including the en- 
tropy term, seem totally missing. This is particularly frustrating in view of 
the desire to characterize and possibly predict surface phase transitions as 
a function of parameters, such as temperature, coverage, external stresses, 
crystal bending, etc. In systems where interatomic forces are well under- 
stood, an alternative methods to predict such phase transition is simulation, 
particularly Molecular Dynamics (MD) simulation. However, besides being 
somewhat less fundamental, the simulation approach usually suffers from 
practical problems, such as size and time limitations, which greatly restrict 
the variety of transitions that can be directly described. 

In the present paper we present a route to calculate directly the surface 
free energy change with respect to one specific external parameter, the cur- 
vature of the underlying solid. The method is again based on Molecular 
Dynamics (MD) simulations, and uses the framework for variable curvature 
MD developed by Passerone et al.^. The surface phases to be compared 
can be either simulated separately, or realized on the two opposite faces of 
the same simulation slab. Since free energies are calculated separately, the 
actual phase transformation is not required to take place spontaneously in 
the simulation - in fact, it must be avoided. Curvature of a crystal plate is a 
standard tool for observing surface stress difference of the two opposite faces 
of the plate [|, 0]. Conversely, inducing curvature of the plate by an external 
bending force will cause the two opposite faces to depart from their original 
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state, and their free energies will generally evolve in opposite directions. 

We have chosen metallic surfaces as our test case, for a variety of rea- 
sons. The surfaces of metals are readily accessible to a variety of surface 
techniques. Especially noble and some transition metals show a variety of 
surface reconstructions where the surface stress plays an important rolep[. 
If a metal slab is bent the surfaces are strained and work is exerted onto 
them. It is plausible that bending could eventually drive phase transitions, 
including changes or removal of the reconstruction. Preliminary results on 
reconstructed Au(lll) indicate that the surface reconstruction, or at least 
some features of it, can in fact be removed by curvature]^ 

2 Method 

Although MD is a very useful technique to study the time evolution of ener- 
gies, temperature, pressure, stress, ... of a system with well defined forces, 
other thermodynamic quantities, like the free energy and entropy, are not 
obtained. The reason is that such quantities are not an average of some 
mechanical entity over the ensemble of the states: they contain information 
on the whole ensemble of states which cannot be directly extracted from the 
time evolution of a sample [^]. Nonetheless in some specific cases the free 
energy variation along a reversible, isothermal path can be determined from 
mechanical quantities by integrating the work done onto the system. 

In the particular case of the slab, where curvature is forced externally 
through bending, the work is given by the integral of the stress tensor with 
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respect to the strain and the volume of the slab. This integral contains 
contributions both from the surfaces and from the bulk of the slab. For our 
purposes, the former must be separated from the latter 0]. The first step is 
to extract the total stress field cTq/j in the simulation of a curved slab. The 
only relevant component is ayy, where the y-axis is parallel to the surface and 
is oriented in the direction that is stretched during bending (the other axes 
are chosen with z normal to the slab). The average force exchanged across 
the a;|/-plane is0 

where the indexes i and j runs over the particles, U is the potential energy as a 
function of the interatomic distances r^j = l^j— | and L^, Ly, indicate the 
lengths of the slab along the three axes (the latter being the slab thickness) . 
Note that this formula does not require pairwise potentials as it might seem. 
It applies for arbitrary many body potentials (including EAM, glue models, 
etc) such that total energy depends only on pair distances. In the following 
for convenience we shall write Fij in place of dU/dvij. Note that this is 
generally not the force acting between particles i and j, but it becomes that 
in the specific case of pairwise additive interactions. The first term in the 
RHS is the kinetic energy contribution to stress, which in the classical case 
is just [NksT /V) ■ L^Ly, being the number of atoms and V = L^LyL^ the 
volume of the cell. The second term arises from interparticle interactions: it 
is the derivative of the potential energy with respect to a uniform stretching 
along the y-axis, averaged over the system configurations. 



We cut a slab in slices along the z direction, each slice corresponding to 
one layer (the depth coordinate z is orthogonal to the bending axis x, and 
to the bending direction y). By restricting the sum to particles located in 
each slice we calculate the stress along the bending direction, resolved layer 
by layer. For a bent slab, equation (|l]) must be generalized to compute the 
average stress per layer in the direction which locally follows the profile of 
the surface. This is done by computing the potential energy derivative for 
a uniform strain of the slab along y. In order to better deal with the sym- 
metry of the bent slab, we introduce a system of scaled, curved coordinates 
(gi,g2,Q'3); indicating with k the slab curvature in the neutral cylinder (i.e. 
in the middle of the slab): 

X = L^qi 

y = {l/k + L^q-i) sm{kLyq2) (2) 
z = {l/k + L^qs) cos{kLyq2) -1/k 
The variable q2 is proportional to the polar angle 6 through 6 = kLyq2] thus q2 
moves along the bending direction. The "radial" coordinate gs is measured 
along the direction normal to the surface and is zero in the middle layer. 
Each of these coordinates ranges between —1/2 and 1/2. In the limit of zero 
curvature they are proportional to x, y and z respectively and correspond 
to the scaled coordinates introduced by Andersen |Tl[] and by Parrinello and 



Rahman ||12||. We note that this system of coordinates is slightly different 
from the one introduce in ref. 0: the present choice is more convenient in 
order to compute the stress field. Omitting the details of the calculation, the 
final formula for the derivative of the potential energy keeping all the scaled 
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coordinates fixed is 



dU 



dL 



Y 
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Since periodic boundary conditions are used, both the distances rij and the 
differences 6ij = 6i~6j = kLy{q2i — q2j) must be computed with the minimum 
image convention. The contribution of particle i to the force J j dxdz 022 is 
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The average force any layer exerts in the direction is obtained as: 
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In the zero curvature hmit , k —>■ 0, this formula reduces to (|l|). 

To exemplify, we apply this scheme to the Au(OOl) surface, modeled by 
means of the "glue" many body potential ||13|. It is well known that this 
surface reconstructs with a denser triangular top layer |14, 1^, increasing its 
lateral density by about 24% relative to an unreconstructed layer. One could 
expect that the tensile stress typical of the unreconstructed surface should 
have decreased, maybe even disappeared, or reversed. Moreover, different 
reconstruction periodicities with different lateral densities and different sur- 
face stresses might come in competition with one another for the lowest free 
energy, as curvature is cranked up. Figure |l] shows the stress T22/LX in 
each slice of a (001) gold slab. The two slab surfaces are given two slightly 
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Figure 1: The linear stress (force/length) associated to every layer of two 
samples with different thicknesses. In the insets the reconstructions of the two 
surfaces are shown. Surface atoms with the largest \z\ are shown as lighter, 
and identify the reconstruction solitons. Note the bulk linear z dependence, 
and the large surface contribution. 
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different reconstructions, namely 1x5 (close to the experimental one) and 
1x6 (with slightly lower surface density). Two slab samples are considered: 
one made by 12 layers and a thicker one made by 24 layers. Both samples 
have sizes L^. = 28. SA and Ly = 172. 9A. The bending direction is [110]. The 
stress is computed through (|^) taking the averages during a 36 psec evolution 
time at T = 100 Kelvin. The bending curvature in the middle plane z = 
where strain vanishes, is k = 4.49 ■ 10~^A ^ . The strain of such a curved 
slab is linear with z by construction and is determined by the curvature k 
through the simple relation £22 = kz. The z-resolved stress distribution is 
instructive. The stress profile deep in the slab is, as expected, again linear 
with z, and proportional to the strain. Close to the surfaces, however, the 
stress differs from its bulk-like extrapolation, and oscillates. Both surface 
top layers show a positive (tensile) stress, that is the surfaces further tend 
to reduce their area. This is rather the rule on metal surfaces]^], but it is 
interesting to note that even the more close-packed 1x5 reconstruction has 
not quite eliminated the tensile stress. The oscillations, in turn, reflect the 
composite layer-dependent response to the main tensile force exerted by the 
top layer. In order to calculate the total surface stress, we must subtract 
from the actual stress distribution the corresponding extrapolated bulk lin- 
ear stress, and integrate. The final result must be independent of the slab 
thickness. Comparison of the slab-resolved stress of the thin slab with that 
of the thicker one, as in Figure |l], shows that this is indeed the case, and 
that the procedure works even in the thinner one. In the following we will 
write T22^(/) for the surface contribution to the force associated to layer I. 
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The surface stress is 

where the sum extends over the layer close to the interface with not negli- 
gibly small ■ The value obtained for Au(lOO) is of 4.0 N/m in excellent 
agreement with the experimental value of 4.4 N/m [|16| . 



3 Curvature-Dependent Free Energy 

Now we are ready for the next step in the free energy difference calculation. 
As the curvature is varied in an isothermal environment the reversible work 
done against the surface contribution should equal the surface free energy 
variation. Since both strain and stress depend on the layer depth (£22(0 — 
kzi), the sum must be done layer by layer 

ciF(^)=L,^T2?(/)&22(0 
I 

Replacing the strain £22(0 with zik and integrating over the curvature k: 

AF(^) = Ly I dk ^T!fXl)zi 
I 

In conclusion the surface free energy per unit area is linked to its value f^^ 
for the flat slab through the relation 

f^'' = fi'^ + ^ fdkY,T^Xl)zi (5) 

We continue our exemplification of the method to the (001) surface of 
gold, where we wish to calculate curvature-dependent free energies. We pre- 
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Figure 2: The free energy per unit area for some reconstructed Au(OOl) 
surfaces. In abscissa the strain at the surface due to bending is reported. 
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pared different (001) slab samples, where the surfaces differ by a different re- 
construction periodicity, in turn related to a different density of added atom 
rows perpendicular to the bending direction, that is parallel to the x-axis. In 
particular, in a (1 x n) reconstruction one extra row is added every n exist- 
ing rows. After the usual equilibration steps the surfaces relax and equally 
spaced solitons appear, corresponding to misfit dislocations between first and 
second layer. The slabs are gradually bent and for each curvature the overall 
layer-resolved stress is computed through (p, the averages taken during a 36 
psec evolution time in a canonical simulation at T = lOOK. From the forces 
T22 the surface contributions are extracted, as explained above. We plot the 
results in Fig. ^ where we have for convenience assumed arbitrarily the un- 
known zero- curvature free energy to be equal to the zero-temperature energy. 



calculated previously through MD by F. Ercolessi et al.ljT^- At our working 
temperature of 100 K, this is likely a not unreasonable approximation. 

All the surface free energies decrease with curvature on the concave side, 
where the compressive subsurface strain acts to reduce the need for a tensile 
surface stress. The decrease is more pronounced, as it should be, for the less 
dense higher reconstruction, whose tensile stress is higher. Although there 
are crossings, no other reconstruction crosses the (1 x 5) surface free energy, 
which is therefore predicted to be stable against curvatures leading to surface 
strains up to the percent range. 

Preliminary as these results are, they seem encouraging, although we do 
not yet have an alternative route to check them. It will be interesting in the 
future to carry out tests aimed among other things at separating the internal 
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energy from the entropy part. 

In conclusion, we have presented a calculation of surface free energy varia- 
tion with curvature, realized by direct integration of surface stress, obtained 
through realistic molecular dynamics simulation. A specific application to 
(001) reconstructed gold surfaces shows a good feasibility of the method, 
and foreshadows future applications to study surface phase transitions under 
curvature. 
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